Integrated bioinformatics analysis of retinal ischemia/reperfusion injury in rats with potential key genes

The tissue damage caused by transient ischemic injury is an essential component of the pathogenesis of retinal ischemia, which mainly hinges on the degree and duration of interruption of the blood supply and the subsequent damage caused by tissue reperfusion. Some research indicated that the retinal injury induced by ischemia-reperfusion (I/R) was related to reperfusion time. In this study, we screened the differentially expressed circRNAs, lncRNAs, and mRNAs between the control and model group and at different reperfusion time (24h, 72h, and 7d) with the aid of whole transcriptome sequencing technology, and the trend changes in time-varying mRNA, lncRNA, circRNA were obtained by chronological analysis. Then, candidate circRNAs, lncRNAs, and mRNAs were obtained as the intersection of differentially expression genes and trend change genes. Importance scores of the genes selected the key genes whose expression changed with the increase of reperfusion time. Also, the characteristic differentially expressed genes specific to the reperfusion time were analyzed, key genes specific to reperfusion time were selected to show the change in biological process with the increase of reperfusion time. As a result, 316 candidate mRNAs, 137 candidate lncRNAs, and 31 candidate circRNAs were obtained by the intersection of differentially expressed mRNAs, lncRNAs, and circRNAs with trend mRNAs, trend lncRNAs and trend circRNAs, 5 key genes (Cd74, RT1-Da, RT1-CE5, RT1-Bb, RT1-DOa) were selected by importance scores of the genes. The result of GSEA showed that key genes were found to play vital roles in antigen processing and presentation, regulation of the actin cytoskeleton, and the ribosome. A network included 4 key genes (Cd74, RT1-Da, RT1-Bb, RT1-DOa), 34 miRNAs and 48 lncRNAs, and 81 regulatory relationship axes, and a network included 4 key genes (Cd74, RT1-Da, RT1-Bb, RT1-DOa), 9 miRNAs and 3 circRNAs (circRNA_10572, circRNA_03219, circRNA_11359) and 12 regulatory relationship axes were constructed, the subcellular location, transcription factors, signaling network, targeted drugs and relationship to eye diseases of key genes were predicted. 1370 characteristic differentially expressed mRNAs (spec_24h mRNA), 558 characteristic differentially expressed mRNAs (spec_72h mRNA), and 92 characteristic differentially expressed mRNAs (spec_7d mRNA) were found, and their key genes and regulation networks were analyzed. In summary, we screened the differentially expressed circRNAs, lncRNAs, and mRNAs between the control and model groups and at different reperfusion time (24h, 72h, and 7d). 5 key genes, Cd74, RT1-Da, RT1-CE5, RT1-Bb, RT1-DOa, were selected. Key genes specific to reperfusion time were selected to show the change in biological process with the increased reperfusion time. These results provided theoretical support and a reference basis for the clinical treatment. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-024-10288-0.


Introduction
As one of the highest oxygen-consuming tissues in the body, the retina has a delicate and complex structure and a vigorous metabolism [1].The retina is a transparent membrane that is divided into ten layers from the outside, namely, the retinal pigment epithelium, the cone and rod layer of the retinal neurosensory layer, the outer limb, the outer nuclear layer, the outer plexiform layer, the inner nuclear layer, the inner plexiform layer, the ganglion cell layer composed of ganglion cell nuclei, the nerve fiber layer composed of nerve cell axons, and the inner limb [2,3].Retinal tissue has a dual blood supply system for adequate energy supply and oxygen consumption: the outer five layers of the retina are supplied by the choroidal blood vessels (referred to as the choroidal circulation) [4][5][6], whereas the inner five layers are supplied by the central retinal artery (referred to as the retinal circulation) [7,8].The choroidal circulation is much denser than the retinal circulation, and the sparse retinal circulation is more favorable for light passage but makes the retina more susceptible to vascular disease [9].Any pathological damage or retinal ischemia and hypoxia caused by retinal vascular obstruction can lead to infarction of the retinal tissue cells, which results in the loss of the function of receiving and transmitting light stimuli from the external environment.After the blood supply is restored, the damage to the structure and function of the retina is further aggravated, resulting in retinal ischemia /reperfusion (RIR) injury [10,11].A series of scientific studies have confirmed that the pathogenesis of RIR injury also includes a series of deleterious events, including complement system activation and leukocyte recruitment, endoplasmic reticulum stress, calcium overload, decreased oxidative phosphorylation, increased free radicals, vascular endothelial cell dysfunction, apoptosis signaling, necrosis and autophagy [12][13][14].It is the main cause and pathogenesis of retinal thinning and atrophy, retinal ganglion cell (RGC) death, and visual function impairment [15].
Non-coding RNA (ncRNA) regulation is one of the essential epigenetic regulatory mechanisms.ncRNAs are functional RNA molecules that cannot be translated into proteins and have regulatory roles, mainly including microRNA (miRNA), long non-coding RNA (lncRNA), and circular RNA (circRNA) [16,17].In recent years, competing endogenous RNA (ceRNA) networks consisting of ncRNA and mRNA have gradually attracted the attention of researchers.The nature of the ceRNA network is that lncRNA or circRNA competes with corresponding mRNA and binds to the same miRNA, thereby participating in the regulation of gene expression at the post-transcriptional level [18,19].The primary ceRNA networks are circRNA-/lncRNA-miRNA-mRNA networks.Findings suggest that ceRNA networks may have critical regulatory roles in the development of RIR processes in multiple organs, including RIR [20][21][22].For example, lncRNA Ttc3-209 is significantly up-regulated after RIR injury, which up-regulates the translation of Wnt8a mRNA through sponging miR-484, thereby promoting RGC cell apoptosis [20].However, screening critical genes and their specific mechanism and functions involved in the development of RIR injury still needs more investigation.
The primary purpose of this study was to screen and identify the key genes in RIR injury by integrating mRNA, lncRNA and circRNA sequencing data, to provide theoretical support and reference basis for clinical treatment.In this study, we screened the differentially expressed circRNAs, lncRNAs and mRNAs between the control group and the model group, and within different reperfusion times (24h, 72h, 7d) using whole-transcriptome sequencing, and obtained the expression trends of the time-varying mRNAs, lncRNAs, and circRNAs by time-sequencing analysis.Then, candidate circRNAs, lncRNAs, and mRNAs were obtained by intersection of differentially expressed genes and tendency change genes.key genes whose expression changed with the prolongation of reperfusion time were selected by the importance scores of the genes.The biological pathways and potential regulatory mechanisms of the key genes were analyzed by bioinformatics, and the drugs associated with them as well as the specific molecular regulatory mechanisms during different periods of reperfusion were investigated.In addition, characteristic differentially expressed genes specific to the reperfusion time were analyzed, and key genes specific to the reperfusion time were screened to show the changes in biological processes with the prolongation of reperfusion time.The workflow of the study was listed in Fig. 1.

Construction of high intraocular pressure (IOP) model
High IOP model was constructed for whole transcriptome sequencing.12 Sprague Dawley (SD) rats (Half male and half female, 220±20 g) were purchased from Beijing Silaikedake Laboratory Animal Technology Co., LTD.Rats were fed and managed in strict compliance with the Vision and Ophthalmology Research Association Statement and approved by the Ethics Committee of the First Affiliated Hospital of Kunming Medical University.SD rats were randomly divided into a control (rats was maintained at normal IPO levels), model_24h group (subjected to high IOP process and followed with reperfusion for 24h), model_72h group (subjected to high IOP process and followed with reperfusion for 72h), and model_7d group (subjected to high IOP process and followed with reperfusion for 7 days), with three rats in each group.As mentioned before [23,24], all groups of rats were injected intraperitoneally with 100 mg/kg ketamine and 5 mg/kg xylazine to anesthetize rats for experimental injury, and 0.5% Alcaine ophthalmic solution was used to anesthetize the corneas of rats topically, and 1% tropicamide was used to dilate the pupils.Subsequently, rats in the model_24 h, model_72 h, and model_7 d groups were subjected to corneal cannulation to elevate IOP to 110 mmHg for 60 min to induce the ischemic process.Corneal cannulation in the control group of rats was maintained at normal IOP levels.Further, the eyes of rats in the model_24 h, model_72 h, and model_7 d groups performed reperfusion for 24 h, 72 h, and 7 d, respectively.Finally, the retinas from each group were collected for whole transcriptome sequencing.

RNA extraction and quality assessment
Total RNA was extracted from the retinal tissues of each group of rats using the mirVana miRNA Isolation Kit (Ambion, 1561) with reference to the manufacturer's Fig. 1 The Workflow of this study protocol.Subsequently, total RNA integrity was qualified and assessed by NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, MA, USA) and Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA).Samples with RNA integrity number ≥7 were subjected to subsequent analysis.

Library construction and high-throughput sequencing
RNA library construction is the process of converting RNA into double-stranded DNA recognizable by a second-generation sequencer through processes such as reverse transcription and splice junctions.Construction of RNA libraries from qualified RNA samples was complicated using TruSeq Stranded TotalRNA with Ribo-Zero Gold (Illumina, CA, USA).For lncRNA sequencing, rRNA was removed by Ribo-Zero Gold, and then libraries of lncRNA were prepared.For circRNA sequencing, rRNA was removed, and RNA was treated with RNase R (Epicentre, Madison, WI, USA) to remove linear transcripts and was then fragmented to approximately 200 bp.The purified RNA fragments were subjected to first-and second-strand cDNA synthesis following adaptor ligation according to the NEB Next Ultra ™ RNA Library Prep Kit manual for Illumina (NEB, USA).Subsequently, these libraries were sequenced using an Illumina Hiseq 2500 sequencing platform (Illumina), and 150 bp paired-end reads were generated.The raw reads in FastQ format were further quality filtered by Trimmomatic [25] software for adapter removal and filtering.Clean reads were quality inspected by fastqc [26] and then aligned to the human reference genome using hisat2 [27].The sequencing reads of each sample was aligned with the mRNA sequence by transcript4sequences, aligned with the known lncRNA sequences, and lncRNA prediction sequences by bowtie2.Data was imported to eXpress to make a gene quantitative analysis.The FPKM value and counts value were obtained.For circRNAs, we used BWA software [28] to align the sequencing reads of each sample with the reference genome.CIRI2 was used to scan the circRNA-paired chiastic clipping signal [29,30].The expression values were normalized by the reads per million (RPM) algorithm [31], and the number of junction reads counts and fold changes were normalized by DEseq2 in the R script.False discovery rate correction (adj.P<0.05) was used in the analysis.

Differentially expressed gene screening
The objective of differential expression analysis is to find genes that have significant differences in expression levels between samples.These samples can represent different biological states, such as drug-treated versus control, diseased versus healthy individuals, different tissues or different developmental stages.Principal component analysis (PCA) was a widely used data dimensionality reduction algorithm to obtain the overall degree of similarity or difference in gene expression values between samples, which was accomplished with scatterplot3d package (v0. .The differentially expressed cir-cRNAs, lncRNAs, and mRNAs between model_24h group and control group, model_72h group and control group, model_7d group, and control group, and model group (the combination of samples in model_24h group, model_72h group, and model_7d group) and control group were screened using the DEseq2 R package.The threshold value was set as |log2 (fold change, FC) | > 0.5 and P adj.< 0.05 for lncRNAs and mRNAs, |log2 (fold change, FC) | > 0.5 and P < 0.05 for circRNAs.The ggplot2 R package was utilized to draw volcano plots of the DEGs.The top 10 up-regulated genes and the top 10 down-regulated genes are labeled in the volcano plot sorted according to the ploidy of difference log2FC Heat maps were constructed using the heatmap package in the R script.

Identification of Chronological Expression Analysis
In order to identify the circRNAs, lncRNAs, and mRNAs that changed in the same trend with the increase of reperfusion time, the gene expression in each group was clustered using the R package Mfuzz based on fuzzy c-means clustering (FCM), with a desired number of clusters of 10.Mfuzz is a software tool that can be used to analyze gene changes over time.It can also help to understand how genes interact with each other in biological processes and how they affect each other's functions.The gene expression in control (with the reperfusion time as 0h), 24 h reperfusion, 72 h reperfusion, and 7 d reperfusion were set as point1, point2, point3 and point4.All the genes in the sequencing results with the same expression trend in different periods were analyzed by clustering, and then the trend genes were merged to obtain the trend change genes.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment
The GO database provides specialized terminology to define the properties of gene products, which encompasses biological processes (BP), cellular components (CC), and molecular functions (MF) [32].BP represents a process of molecular activity events, including a functional collection of cells, tissues, organs, and species.CC indicates the cell or the external environment it is in.MF describes the active component of the gene product at the molecular level.KEGG is a database that systematically analyzes gene function and links genomic and functional information [33].In this study, GO and KEGG pathway enrichment analysis of DEmRNAs was performed by the clusterProfiler package [34].P <0.05 is considered statistically significant.

Protein-protein interaction network (PPI)
The PPI of candidate mRNAs and time-specific mRNAs were obtained by String (https:// string-db.org/).Cystoscope was utilized to visualize the interconnections between differential genes using the degree of node importance (Degree) as a metric, which indicates the number of neighbors of the node, and highly connected nodes indicated a more critical role in the network.The CytoHubba was used to calculate the importance of the genes in the modules using the topology algorithm MCC.

Gene Set Enrichment Analysis (GSEA) enrichment analysis for candidate genes
GSEA was performed to detect the biological function, chromosomal location, or regulation of genes [35].GSEA is a gene function analysis method that is used to reveal differences in gene expression under two different conditions.GSEA focuses on the co-regulatory patterns of an entire gene set, rather than on the expression levels of individual genes.Through statistical methods, GSEA is able to identify the enhancement or inhibition of a biological pathway or function, thereby revealing the underlying biological mechanism.Here, the correlation between key genes and other genes in the transcriptome sequencing was determined by Pearson correlation analysis to rank genes.The CP: KEGG dataset under C2 in the rat (Rattus norvegicus) species was extracted as a background gene set by the R package msigdbr and analyzed for GSEA enrichment.Finally, GSEA enrichment analysis was performed using the GSVA R package.The enrichplot R package was used to visualize the Top 5 pathways.

Regulatory networks and prediction
The competitive endogenous RNA (ceRNA) regulatory network consists of mRNAs, coding pseudogenes, longchain noncoding RNAs, and miRNAs.ceRNA networks form a highly interconnected regulatory system by mediating the competitive binding of miRNAs between different RNA molecules.The significance of studying ceRNA networks lies not only in expanding the knowledge of gene regulatory mechanisms, but also in their potential applications in various biological and medical fields.To construct ceRNA networks, the miRNAs targeted to key genes were predicted by miRWalk (Ver.2.0, http:// mirwa lk.umm.uni-heide lberg.de/).And the binding relationship of candidate lncRNAs and candidate circRNAs with miRNA was predicted by miRanda (http:// www.micro rna.org/ microrna/home.do).To find other genes associated with the key genes, GeneMANIA (http:// www.genem ania.org/) was used to predict the genes associated with the biomarker function and the functions involved.Subcellular localization pertains to the specific site where a protein or its expression product is found within a cell, such as the nucleus, cytoplasm, or cell membrane.Increasing evidence indicates that RNAs residing in various subcellular organelles exhibit distinct functionalities during biological processes, highlighting the significance of subcellular localization in unraveling the intricate biological functions of RNAs.For subcellular location, the base sequences of five key genes (Cd74, RT1-Da, RT1-CE5, RT1-Bb, RT1-DOa) were queried through the NCBI website (https:// www.ncbi.nlm.nih.gov/).The base sequences of the genes were then entered into the mRNALocater (http:// bio-bigda ta.cn/ mRNAL ocater) database to obtain predicted scores for the five subcellular localizations, with the highest score being used for the final specific, accurate localization.X2K (eXpression2Kinases, https:// amp.pharm.mssm.edu/ X2K/) is used to identify upstream regulators that may be responsible for the observed genome-wide gene expression patterns by integrating information based on transcriptome data and existing signal transduction databases.To explore the potential regulatory mechanisms of key genes, proteins that interact with transcription factors and upstream kinases that regulate key genes were analyzed by eXpres-sion2Kinases.SIGNOR (http:// signor.uniro ma2.it/) is a manually annotated database of causal relationships between human proteins, biologically relevant chemicals, stimuli, and phenotypes.In order to analyze the signaling relationships of critical genes further, key genes were annotated through the SIGNOR database to construct an essential gene signaling information network.The DGIdb (https:// www.dgidb.org) database was utilized to search for drugs targeting key genes for potential drug discovery of biomarkers and to construct a drug-gene interaction network.Since the database is for human-targeted genes, we used the R package "homologene" to convert rat key genes into homologous human genes and then searched the DGIdb website.To analyze the role of key genes in other eye diseases, the Comparative Toxicogenomics Database was used (CTD, https:// ctdba se.org/) to analyze the relationship between ocular diseases and key genes.Bar charts were drawn to show the top five diseases with the highest predicted scores for each key gene respectively.

Real-tine quantitative reverse transcription (RT-qPCR) assay
The mRNA expression of key genes was detected by RT-qPCR assay.The extracted total RNA was used as a template, and reverse transcription was carried out using an RNA first-strand cDNA synthesis kit.The relative expression levels of mRNA were determined by the qPCR method.qPCR reaction conditions were as follows: 95℃ 3 min, 95℃ 15s, 60℃ 35s, 72℃ 30 s, a total of 40 cycles.GAPDH was used as the internal reference of mRNA, and the relative expression level was calculated by the 2 -ΔΔCt method.The primers were listed in Table 1.

Statistical analysis
GraphPad Prism 8.0 was used to analyze the data and graph.Data were presented as mean ± standard deviation ( x ± SD ) and Turkey's test was used for comparison between the two groups.P<0.05 represents the difference is significantly significant.

Chronological analysis of differentially expressed genes
On the other side, chronological changes in gene expression with increasing reperfusion time were detected by the Mfuzz algorithm.The R package Mfuzz was used to perform cluster analysis of the model and control groups to identify the RNAs with trend changes in samples from different reperfusion time.All the RNAs obtained from the sequencing results with the same expression trend in different reperfusion periods (control group, model_24h group, model_72h group, and model_7d group) were clustered.Then, the trend RNAs were concatenated to obtain the trend change RNAs.As a result, mRNA expression in cluster 1 and 7 were increased with the increase of reperfusion time, and therefore, the genes in the clusters were merged, and 3733 trend mRNAs were obtained (Fig. 3A, Supplementary Tables 13 and  14).The lncRNA expression in cluster 5 and 9 were increased, while lncRNA expression in cluster 3 were decreased with the increase of reperfusion time, and  4584 trend lncRNAs were obtained (Fig. 3B, Supplementary Tables 15 and 16).The circRNA expression in cluster 7 and 9 were increased, while circRNA expression in cluster 8 and 10 were decreased with the increase of reperfusion time, and 4635 trend circRNAs were obtained (Fig. 3C, Supplementary Tables 17 and 18).Then, 316 candidate mRNAs, 137 candidate lncRNAs, and 31 candidate circRNAs were obtained by the intersection of differentially expressed mRNAs, lncRNAs, and circRNAs with trend mRNAs, trend lncRNAs and trend circRNAs (Fig. 4A-C, Supplementary Tables 20 and 21).To evaluate the biological functions that candidate genes enriched, GO annotation and KEGG enrichment was preformed.Candidate mRNAs were enriched in 686 GO terms (629 in BP terms, 21 in CC terms, and 36 in MF terms), such as positive regulation of cell activation, adaptive immune response, positive regulation of leukocyte activation, external side of plasma membrane, extracellular matrix, external encapsulating structure, immune receptor activity, cytokine binding, peptide binding, involving genes such as A2m, Acp5, C1qc, Acta2, Bgn, Cd180, Cd1d1, Cd74, Clu, Also, candidate mRNAs were enriched in 57 KEGG terms, such as Cell adhesion molecules, Cytokinecytokine receptor interaction, Epstein-Barr virus infection, Human T-cell leukemia virus 1 infection, Coronavirus disease-COVID-19, including genes such as RT1-Da, Icoslg, RT1-Db1, RT1-Ba, RT1-Bb, Icos (Fig. 4D-F, Table 2, Supplementary Tables 22 and 23).

Key gene identification and the prediction of regulation
In order to understand the interactions between proteins encoded by candidate genes, 316 candidate mRNAs were uploaded to the STRING database for PPI network construction (interaction score >0.7), the protein-protein interaction network that contained 297 nodes and 281 edges was obtained (average node degree: 1.89).Then, PPI was imported to Cytoscape (Ver.3.7.2),and the importance scores of the genes were obtained by the topology algorithm MCC to identify key genes.
Then, GSEA was used to evaluate the function of key genes.Cd74 was found to be positively related to ribosome, antigen processing, and presentation, regulation of actin cytoskeleton (Fig. 6a, Supplementary Table 34).RT1-Bb was found to be positively related to ribosome, antigen processing, and presentation (Fig. 6b, Supplementary Table 26).RT1-CE5 was found to be positively related to ribosome (Fig. 6C, Supplementary Table 27).RT1-Da was found to be positively related to ribosome, antigen processing, and presentation, regulation of actin cytoskeleton (Fig. 6D, Supplementary Table 28).RT1-DOa was found to be negatively related to ribosome (Fig. 6E, Supplementary Table 29).
To predicted the targeted miRNAs and ncRNA with key genes, ceRNA networks was constructed.A network included 4 key genes (Cd74, RT1-Da, RT1-Bb, RT1-DOa), 34 miRNAs and 48 lncRNAs and 81 regulatory relationship axes (32 related with Cd74, 12 related with RT1-Da, 23 related with RT1-Bb, 14 related with RT1-DOa) was constructed (Fig. 7A, Supplementary Table 31).Also, a network included 4 key genes (Cd74, RT1-Da, RT1-Bb, RT1-DOa), 9 miRNAs and 3 circRNAs (circRNA_10572, circRNA_03219, circRNA_11359) and 12 regulatory relationship axes (4 related with Cd74, 1 related with RT1-Da, 5 related with RT1-Bb, 2 related with RT1-DOa) was constructed (Fig. 7B, Supplementary Table 32).RT1-CE5 was not included in the networks since no miRNA targeted was predicted in the miRWalk database.GeneMA-NIA was utilized to predict the genes associated with the key genes and the functions involved to find the genes associated with the key genes.20 related genes with 1144 total links were predicted (Physical Interactions accounted for 77.64%,Co-expression 8.01%, Predicted 5.37%, Co-localization 3.63%, Genetic Interactions 2.87%, Pathway 1.88%.Shared protein domains accounted for 0.6%).These related genes were enriched in 98 functions, including antigen processing and presentation of peptide antigen, antigen processing and presentation, MHC protein complex, antigen binding, antigen processing and presentation of exogenous peptide antigen, antigen processing and presentation of exogenous antigen, antigen processing and presentation of peptide antigen via MHC class II (Fig. 7C, Supplementary Tables 33 and 34).
mRNALocater was used to predict the subcellular location of key genes.Based on the 5 highest scores for subcellular localization, RT1-CE5 localized to the nucleus, Cd74, RT1-Da, RT1-Bb, and RT1-DOa localized to the cytoplasm (Fig. 7D, Supplementary Table 35).In order to explore the potential regulatory mechanisms of key genes, transcription factors, proteins that interact with transcription factors, and upstream kinases that regulate key genes were analyzed by eXpression2Kinases.The results showed that a total of 10 transcription factors (ZNF778, AEBP2, GCM2, ADNP, ESR2, ZNF577, NR3C2, ZNF516, MESP2, HES1), 28 proteins interacting with transcription factors, and 10 upstream kinases (BUB1B, CDK2 CSNK2A1, AKT1, HIPK2.PINK1, AURKB, ABL1, CDK7, PDPK1) were predicted (Fig. 7E, Supplementary Table 36).In order to analyze the signaling relationship of key genes further, the key genes were annotated by the SIGNOR database to construct the key gene signaling network.The results showed that the proteins CIITA and RFX5 were transcriptional activators of RT1-Da, the complex RFX was also transcriptional activator of RT1-Da, the protein NFX1 was transcriptional repressor of RT1-Da, the proteins Marchf1 and March 8 were repressor of RT1-Da, and the complex EBVgH:gL: gp42 was activator of RT1-Da.gp42 activates RT1-Da protein.At the same time, RFX also activates RT1-Bb and RT1-DOa proteins (Fig. 7F, Supplementary Table 37).
Also, DGIdb predicted the drugs that targeted key genes, and the network was constructed.A total of 27 drugs were predicted to be targeted with 4 key genes (Cd74, RT1-Da, RT1-CE5, RT1-Bb), Cd74 was predicted as the target of 1 drug (MILATUZUMAB), RT1-Da was predicted as the target of 6 drugs, RT1-CE5 was predicted as the target of 15 drugs, RT1-Bb was predicted as the target of 12 drugs.Also, 2 drugs (CLAVULANIC ACID, AMOXICILLIN) were predicted to be common targeted to RT1-Da, RT1-CE5, and RT1-Bb. 2 drugs (CARBAMAZEPINE, TICLOPIDINE) were predicted to be common targeted to RT1-CE5, RT1-Bb. 1 drug Fig. 5 The protein-protein interaction of candidate mRNAs.Left: key genes obtained by MCC algorithm; Right: The protein-protein interaction of candidate mRNAs  38).
To analyze the role of key genes in eye diseases, the relationship between eye diseases and key genes was analyzed using the Comparative Toxicogenomics Database, and bar graphs were plotted to show the top five diseases with the highest predicted scores for each key gene, respectively.The results showed that the Cd74 gene was mainly associated with diseases such as Retinal Diseases, Eye Abnormalities, Cataract, and Vision Disorders.RT1-Bb Gene was mainly associated with Cataract, Retinal Degeneration, Eye Abnormalities, and Dry Eye Syndromes.RT1-CE5 genes were mainly associated with Eye Abnormalities, Pathologic Nystagmus, Vision Disorders, Retinoschisis, and Horner Syndrome.The RT1-Da gene was mainly associated with Vision Disorders, Optic Nerve Diseases, Cataract, Retinal Diseases, and Eye Abnormalities.The RT1-Da gene was mainly associated with Vision Disorders, Optic Nerve Diseases, Cataract, Retinal Diseases, and Eye Abnormalities.The RT1-DOa gene was mainly associated with Retinal Diseases, Lens Subluxation, Night Blindness, Uveitis, and Anterior and Corneal Injuries.In conclusion, key genes were mainly associated with Retinal diseases, Eye abnormalities, Cataract, and Vision disorders (Fig. 7H, Supplementary Table 39).In order to observe the changes in the expression of key genes over time, vertical scatter plots were plotted to show the expression levels of the key genes in The expression of the key genes was detected in transcriptome data and verified by qPCR to observe the changes in the expression among different reperfusion time.The results showed that the expression of Cd74, RT1-Da, and RT1-CE5 was all increased after reperfusion for 24h, 72h, and 7d compared to the control group.RT1-Bb was significantly higher expressed in model_72h group and model_7d compared to the control group.RT1-DOa was significantly higher expressed in model_7d compared to the control group (Fig. 8 A-J).

Discussion
Retinal ischemia/reperfusion (RIR) injury is a remarkably complex pathophysiological process that is widely seen in a variety of ocular diseases, such as retinal vascular occlusions, glaucoma, diabetic retinopathy, and retinopathy of prematurity, which can lead to blindness [15,36].The nature of RIR is that the blockage and the subsequent restoration of blood flow to the tissues induces a series of oxidative stress and inflammatory effects, which ultimately leads to damage to retinal nerve cells, especially retinal ganglion cells (RGCs) [11,15,37].Recently, the pathogenesis of RIR has been poorly understood, and no clinically approved drugs can effectively rescue ischemic retinal nerve cells.Currently, treatments for RIR damage include the counteraction of oxygen radicals and oxidative stress, inhibition of calcium overload, inhibition of apoptosis, inhibition of inflammatory responses and reduction of retinal edema, and counteraction of neurotoxicity of nitric oxide and excitatory amino acids [11,[38][39][40].Despite the many methods of treating RIR damage, treatment outcomes remain suboptimal.The tissue damage caused by transient ischemic injury is an essential component of the pathogenesis of retinal ischemia, which mainly hinges on the degree and duration of interruption of the blood supply and the subsequent damage caused by tissue reperfusion [41].Some research indicated that the retinal injury induced by ischemia/reperfusion (I/R) was related to reperfusion time.Wang et al. found that retinal edema was seen in the early stage and followed by retina atrophied gradually in 72 h and 144 h as reperfusion time [42].Zhang et al. found loss of cells in the retinal ganglion cell layer was apparent 2 days after I/R injury, and the number of degenerated capillaries increased greatly by 7 to 8 days after the injury [43].These available studies show that the pathological changes in the retina after different reperfusion times vary considerably.This may indicate that the treatment strategy for retinal ischemia should be adapted according to the reperfusion time.
In this manuscript, we screened the differentially expressed circRNAs, lncRNAs, and mRNAs between the  , 72h, and 7d) with the aid of whole transcriptome sequencing technology, and the trend changes in time-varying mRNA, lncRNA, circRNA were obtained by chronological analysis.Then, candidate circRNAs, lncRNAs, and mRNAs were obtained as the intersection of differentially expression genes and trend change genes.The candidate mRNAs were mainly enriched in immunerelated terms such as lymphocyte-mediated immunity, complement and coagulation cascades, Hematopoietic cell lineage, antigen processing, and expression.The blood-retinal barrier consists of the inner tight junction between retinal capillary endothelial cells and the outer tight junction between the RPE [44,45].On the one hand, due to the structure of the blood-retinal barrier, the macromolecular antibodies in the retinal vessels and choroidal vessels cannot play their functions.On the other hand, there are no lymphoid tissues in the retina, so the antigens do not cause the clonal proliferation of specific T cells or B cells.The retina, therefore, has long been recognized as a privileged site for immunity [11].However, the immunologic response to various stress cues has been found to play a pivotal role in the retina [13].When ischemia and reperfusion occur, the permeability of the blood-retinal barrier is changed, and the activation of microgila [46,47] and the increase of Treg cells [48,49] are observed.Toll signaling activation was also found and induced inflammasome formation [50].Our results of candidate gene enrichment also showed that immune response was crucial in the process of RIR.
Importance scores of the genes selected the key genes whose expression changed with the increase of reperfusion time.As a result, 5 key genes, Cd74, RT1-Da, RT1-CE5, RT1-Bb, and RT1-DOa, were selected.Cd74 is a receptor for the cytokine macrophage migration inhibitory factor (MIF) [51].CD74 regulates T-cell and B-cell development, dendritic cell (DC) motility, macrophage inflammation, and thymic selection [52].Cd74, RT1-Ba, RT1-Bb, RT1-Da, and RT1-Db1 were referred to as major histocompatibility complex (MHC) class II members [53].The MHC-II molecule is a central molecule in the protein presentation pathway.It binds to processed short peptides and presents them to T lymphocytes, activating them to become effector T cells [54].Abcouwer et al. [55] Minocycline was particularly effective in decreasing the appearance of MHCII+ inflammatory leukocytes and reduced leukocyte adhesion and invasion, as well as vascular permeability in RIR.The result of GSEA showed that key genes were found to play vital roles in antigen processing and presentation, regulation of the actin cytoskeleton, and the ribosome.The research from Honjo [56] showed that the cytokines and growth factor in the aqueous humor activate Rho after the ischemia happened, and the Rho/ROCK signal transduction participates in RIR injury via rearrangement of the actin cytoskeleton that was attributed to improved outflow.With the increase in reperfusion time, key gene expression also increased.Also, key genes were mainly associated with Retinal diseases, Eye abnormalities, Cataract, and Vision disorders.Clavulanic acid and Amoxicillin were predicted to be commonly targeted to RT1-Da, RT1-CE5, and RT1-Bb, which further proved that the key genes might participate in the process of RIR injury by involving the adaptive immune response.
Also, the characteristic differentially expressed genes specific to the reperfusion time were analyzed.Key genes specific to reperfusion time were selected to show the change in biological process with the increased reperfusion time.Intriguingly, the specific genes in 24h as reperfusion time were mainly enriched in the pathways that related to retinal function and cellular response to external stress, such as camera-type eye development, retina development in camera-type eye, visual perception, detection of abiotic stimulus, detection of light stimulus.However, the specific genes in 72h as reperfusion time were enriched in the pathways related to membrane potential and neurodevelopment, such as regulation of postsynaptic membrane potential, neurotransmitter receptor activity, and ion channel activity.The specific genes in 7 days as reperfusion time were enriched in the pathways related to immune responses, such as adaptive immune response, lymphocyte differentiation, and positive T cell selection.The retina is a highly specialized neural tissue that continues the central nervous system.The existing research on the novel drug for RIR injury was mainly targeted to the retinal ganglion cell damage and inflammation.Guan et al. [57] Puerarin can ameliorate RIR injury by suppressing apoptosis and TLR4/ NLRP3 inflammasome activation in RGCs. Lee et al. [58] proved Nicotinamide mononucleotide significantly suppressed retinal functional damage, as well as inflammation.Our results indicated that the biological process in different reperfusion time seems to be different.This may help reveal the mechanisms of the onset and progression of RIR injury and offer a novel aspect of its treatment.
In summary, we screened the differentially expressed circRNAs, lncRNAs, and mRNAs between the control and model groups and at different reperfusion time (24h, 72h, and 7d). 5 key genes, Cd74, RT1-Da, RT1-CE5, RT1-Bb, RT1-DOa, were selected.Key genes specific to reperfusion time were selected to show the change in biological process with the increased reperfusion time.These results provided theoretical support and a reference basis for the clinical treatment.experimental proposal and did not participate in this manuscript's writing, review, and publication.No economic interests were appeared.

Fig. 2
Fig. 2 Differentially expressed mRNAs, lncRNAs, circRNA in high IOP model.A The volcano plot revealed DEmRNAs between model_24 h and control group.B The heatmap revealed DEmRNAs between model_24 h and control group.C The volcano plot revealed DElncRNAs between model_24 h and control group.D The heatmap revealed DElncRNAs between model_24 h and control group.E The volcano plot revealed DEcircRNAs between model_24 h and control group.F The heatmap revealed DEcircRNAs between model_24 h and control group.G The volcano plot revealed DEmRNAs between model_72 h and control group.H The heatmap revealed DEmRNAs between model_72 h and control group.I The volcano plot revealed DElncRNAs between model_72 h and control group.J The heatmap revealed DElncRNAs between model_72 h and control group.K The volcano plot revealed DEcircRNAs between model_72 h and control group.L The heatmap revealed DEcircRNAs between model_72 h and control group.M The volcano plot revealed DEmRNAs between model_ 7d and control group.N The heatmap revealed DEmRNAs between model_ 7d and control group.O The volcano plot revealed DElncRNAs between model_ 7d and control group.P The heatmap revealed DElncRNAs between model_ 7d and control group.Q The volcano plot revealed DEcircRNAs between model_ 7d and control group.R The heatmap revealed DEcircRNAs between model_ 7d and control group

Fig. 4
Fig. 4 Candidate gene selection and their function enrichment.A Candidate mRNAs were obtained by the intersection of differentially expressed mRNAs, with trend mRNAs.B Candidate lncRNAs were obtained by the intersection of differentially expressed lncRNAs, with trend lncRNAs.C Candidate circRNAs were obtained by the intersection of differentially expressed circRNAs, with trend circRNAs.D Scatterplot showed GO enrichment results of candidate mRNAs.Horizontal coordinate is the enrichment factor, vertical coordinate is the name of the enriched pathway, dot size indicates the number of differential genes enriched into the pathway, color indicates the range of p.adjust.E Chord diagrams of KEGG terms of candidate mRNAs.The color of the left gene ribbon represents the logFC of the gene, and different ribbons on the right represent different pathways.F scatterplot showed KEGG enrichment results of candidate mRNAs.Horizontal coordinate is the enrichment factor, vertical coordinate is the name of the enriched pathway, dot size indicates the number of differential genes enriched into the pathway, color indicates the range of p. adjust

Fig. 7
Fig. 7 Regulation prediction of key genes.A Key mRNA-miRNA -lncRNA Networks.Red dots are mRNAs, yellow triangles are miRNAs, and blue squares are lncRNAs.B Key mRNA-miRNA-circRNA network.Red dots are mRNAs, yellow triangles are miRNAs, and green squares are circRNAs.C GeneMANIA Network.The middle is the key gene, the outer circle is the related genes with similar functions to the key gene, different colors connecting lines indicate different networks, and different colors in the pie chart indicate different functions of the genes.D Subcellular localization of key genes.Horizontal coordinates are genes, vertical coordinates are scoring at different sites, with higher scores indicating a higher likelihood of being at that site.E Regulatory network of transcription factors, kinases and related proteins.Transcription factors in red, kinases in green, and proteins in gray.F SIGNOR-based signaling information network for key genes.G Gene-Drug Network.Genes are shown in red and drug names in blue.H Gene-Disease Prediction.Horizontal coordinate is the prediction score, vertical coordinate is the predicted disease name, and different color bars indicate different key genes

Fig. 8 Fig. 9
Fig.8Expression of key genes.A-E Expression of key genes across groups in transcriptome data.F-J Expression of key genes across groups were detected in the retinal samples by qPCR

Fig. 10
Fig. 10 Function enrichment of time-specific differentially expressed gene.A Scatterplot showed GO enrichment results of 24h_specific mRNAs.Horizontal coordinate is the enrichment factor, vertical coordinate is the name of the enriched pathway, dot size indicates the number of differential genes enriched into the pathway, color indicates the range of p.adjust.B Chord diagrams of KEGG terms of 24h_specific mRNAs.The color of the left gene ribbon represents the logFC of the gene, and different ribbons on the right represent different pathways.C scatterplot showed KEGG enrichment results of 24h_specific mRNAs.Horizontal coordinate is the enrichment factor, vertical coordinate is the name of the enriched pathway, dot size indicates the number of differential genes enriched into the pathway, color indicates the range of p.adjust.D Scatterplot showed GO enrichment results of 72h_specific mRNAs.Horizontal coordinate is the enrichment factor, vertical coordinate is the name of the enriched pathway, dot size indicates the number of differential genes enriched into the pathway, color indicates the range of p.adjust.E Chord diagrams of KEGG terms of 72h_specific mRNAs.The color of the left gene ribbon represents the logFC of the gene, and different ribbons on the right represent different pathways.F scatterplot showed KEGG enrichment results of 72h_specific mRNAs.Horizontal coordinate is the enrichment factor, vertical coordinate is the name of the enriched pathway, dot size indicates the number of differential genes enriched into the pathway, color indicates the range of p.adjust.G Scatterplot showed GO enrichment results of 7d_specific mRNAs.Horizontal coordinate is the enrichment factor, vertical coordinate is the name of the enriched pathway, dot size indicates the number of differential genes enriched into the pathway, color indicates the range of p.adjust.H Chord diagrams of KEGG terms of 7d_specific mRNAs.The color of the left gene ribbon represents the logFC of the gene, and different ribbons on the right represent different pathways.I scatterplot showed KEGG enrichment results of 7d_specific mRNAs.Horizontal coordinate is the enrichment factor, vertical coordinate is the name of the enriched pathway, dot size indicates the number of differential genes enriched into the pathway, color indicates the range of p.adjust

Fig. 11
Fig. 11 Analysis of the regulatory mechanism between the model_24h group and the control group of the model.A The protein-protein interaction of 24h_specific mRNAs.B key genes obtained by MCC algorithm.C Key mRNA-miRNA -lncRNA Networks.Red dots are mRNAs, yellow triangles are miRNAs, and green squares are lncRNAs.D Key mRNA-miRNA-circRNA network.Red dots are mRNAs, yellow triangles are miRNAs, and blue squares are circRNAs

Fig. 12
Fig. 12 Analysis of the regulatory mechanism between the model_72h group and the control group of the model.A The protein-protein interaction of 72h_specific mRNAs.B key genes obtained by MCC algorithm.C Key mRNA-miRNA -lncRNA Networks.Red dots are mRNAs, yellow triangles are miRNAs, and green squares are lncRNAs.D Key mRNA-miRNA-circRNA network.Red dots are mRNAs, yellow triangles are miRNAs, and blue squares are circRNAs

Fig. 13
Fig. 13 Analysis of the regulatory mechanism between the model_7d group and the control group of the model.A The protein-protein interaction of 7d_specific mRNAs.B key genes obtained by MCC algorithm.C Key mRNA-miRNA -lncRNA Networks.Red dots are mRNAs, yellow triangles are miRNAs, and green squares are lncRNAs.D Key mRNA-miRNA-circRNA network.Red dots are mRNAs, yellow triangles are miRNAs, and blue squares are circRNAs

Table 1
Primer sequences for qPCR assay

Table 2
Top 10 GO enrichment and Top 5 KEGG enrichment of Candidate mRNAs

Table 3
Top 10 GO enrichment and Top 5 KEGG enrichment of characteristic differentially expressed mRNA in 24h reperfusion (spec_24h mRNA)

Table 4
Top 10 GO enrichment and Top 5 KEGG enrichment of characteristic differentially expressed mRNA in 72h reperfusion (spec_72h mRNAs)

Table 5
Top 10 GO enrichment and Top 5 KEGG enrichment of characteristic differentially expressed mRNA in 7d reperfusion (spec_7d mRNAs)